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General amplitude equations for reaction-diffusion systems 
near to the soft onset of birhythmicity described by a super- 
critical pitchfork-Hopf bifurcation are derived. Using these 
equations and applying singular perturbation theory, we show 
that stable autonomous pacemakers represent a generic kind 
of spatiotemporal patterns in such systems. This is verified 
by numerical simulations, which also show the existence of 
breathing and swinging pacemaker solutions. The drift of 
self-organized pacemakers in media with spatial parameter 
gradients is analytically and numerically investigated. 
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Oscillatory reaction-diffusion systems, such as the 
Belousov-Zhabotinsky (BZ) chemical reaction, exhibit a 
rich variety of nonlinear wave patterns. The first complex 
pattern discovered in this reaction was the target pattern, 
where concentric waves were emitted by a pacemaker rep- 
resenting a periodic wave source jjj . Subsequently, sim- 
ilar target structures have been observed in many other 
chemical, physical, and biological systems [0-§. A sim- 
ple theoretical explanation of target patterns in oscilla- 
tory chemical systems is that their pacemakers are cre- 
ated by impurities which increase the local oscillation 
frequency in the medium f|. Most of the target pat- 
terns seen in the BZ reaction are indeed caused by small 
local inhomogeneities, such as dust particles. A general 
question is whether self-organized target patterns, rep- 
resenting an intrinsic dynamical property, are also pos- 
sible in reaction-diffusion systems. Examples of stable 
autonomous pacemakers with localized or extended wave 
patterns are indeed known for several reaction-diffusion 
models @-^]. 

The aim of the present Letter is to show that sta- 
ble autonomous pacemakers with extended wave patterns 
represent a generic pattern-forming object in oscilla- 
tory reaction-diffusion systems near the onset of birhyth- 
micity. Birhythmicity, the coexistence of two stable limit 
cycles corresponding to uniform oscillations with differ- 
ent frequencies, is possible in various systems fi~7|-|20|| , in- 
cluding glycolytic oscillations pi]] and the photosensitive 
BZ reaction p3). Here, we derive two coupled amplitude 
equations yielding the normal form of such a dynami- 
cal system near a supercritical pitchfork-Hopf bifurca- 
tion which leads to birhythmicity. Using singular per- 
turbation theory, an analytical solution for autonomous 
pacemakers is then constructed and its stability is nu- 
merically confirmed. In addition, target patterns with 



breathing or swinging pacemakers arc observed. Finally, 
we show that autonomous pacemakers can drift under 
the influence of a parameter gradient and determine the 
drift velocity. 

The derivation of normal forms for various kinds of 
reaction-diffusion systems has recently been discussed by 
one of the authors 23 24]. Here, we focus our atten- 
tion on systems near to a supercritical pitchfork-Hopf 
bifurcation, where a real uniform eigenmode and a pair 
of complex conjugate uniform eigenmodes start to grow 
simultaneously. This implies that either a stable small- 
amplitude limit cycle becomes unstable and gives rise to 
two stable limit cycles or a pair of stable and unstable 
limit cycles of small amplitude emerges near to a stable 
limit cycle. At least three species are needed to realize 
this bifurcation. 

Using the approach described in p3[ , we have derived 
the normal form of this bifurcation for a general reaction- 
diffusion system p5fl . The normal form is given by 

d t A = A- (1 + ia)\A\ 2 A + (1 + i(3)V 2 A + (1 -ie)Az, 

(la) 

(lb) 



T d t z = a - -f\A\ 2 + z - vz A + X 2 V 2 z, 



The system (|l]) represents a complex Ginzburg-Landau 
equation (CGLE) for a supercritical Hopf bifurcation (|Ta|) 
which is coupled to an equation describing an imperfect 
pitchfork bifurcation (lb). Here A is the complex oscil- 
lation amplitude and z is the amplitude of the slow real 
mode. The coefficients r and A, respectively, are the ra- 
tios of the characteristic time and length scales of the 
real and the oscillatory mode. The parameter e specifies 
the frequency shift of the oscillatory mode due to cou- 
pling to the real mode, 7 characterizes the strength of 
the feedback from the oscillatory to the real mode, and 
v determines the nonlinear saturation of the real mode. 



When cr = 0, Eq. (lb) describes a supercritical pitchfork 
bifurcation, whereas a ^ corresponds to an imperfect 
pitchfork (or "cusp" |2g]) bifurcation. Only positive pa- 
rameters 7, v, and a will be considered in this Letter. 
We assume that uniform oscillations are modulationally 
stable in this system, i.e. the Benjamin- Feir-Newell con- 
dition 1 + a/3 > is satisfied, and that the waves have 
positive dispersion, i.e. (3 — a > 0. 

It can be seen that Eqs. ([[]) have solutions corre- 
sponding to two uniform stable limit-cycle oscillations 
with frequencies Sli and and to an unstable limit 
cycle with frequency 0,%. The frequencies are given by 
^1,2,3 = a + (a + e) 21,2,3 where 21,2,3 are the real roots 
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of the equation vz z — (1 — 7)2 + 7 = a. When a + e < 
(which is the case that we chose in the simulations), the 
smallest root z\ corresponds to the most rapid oscilla- 
tions, i.e. f?3 < fi 2 < 

In addition to these uniform oscillations, the system 
described by Eqs. (0) may have stable nonuniform solu- 
tions representing self-organized pacemakers. To create a 
pacemaker, a sufficiently strong local perturbation should 
be applied to the state comprised of uniform oscillations 
with the lower frequency Q3, so that a small core region 
is formed where oscillations have a higher frequency. In- 
side this region, the variable z is close to z\, whereas 
outside, it is near to z 3 . The core starts to send out 
waves and hence a pacemaker is created. The core ex- 
pands (7 > a is a necessary condition for this, otherwise 
it will contract) and the frequency and the wavenumbcr 
of the emitted waves slowly increase at the same time. In 
turn, this leads to a decrease of the oscillation amplitude 
of emitted waves. This amplitude controls the propa- 
gation velocity of the front, representing the boundary 
of the expanding core. When a critical wavenumber is 
reached, the front velocity becomes zero and a station- 
ary pacemaker is formed. 

We have constructed an analytical solution for sta- 
tionary pacemakers in the one-dimensional system (0). 
The phase <f> and amplitude p are introduced by 
A = pexpj— i(Clst + (/))], Then p is adiabatically elim- 
inated using p 2 ss 1 + z - (V0) 2 + /3V 2 </>. The phase 
equation approximation p?j | is valid for smooth phase 
perturbations and v ;» 1 — 7. Assuming that the char- 
acteristic length scale of the real mode, determining 
the front width, is much shorter than the characteristic 
length scale of the oscillatory subsystem (i.e. that A <C 1), 
we apply singular perturbation theory to this problem. 
The derivation will be published separately and only 
selected results are reported in this Letter. 

The velocity of the front V (of the expanding core) 
depends on the wavenumber fc and is given by 



V(k) = 3~J^z 2 (k), 



(2) 



where z% (k) is the middle root of the cubic equa- 



tion vz 6 — (1 — 7 + 7a) z + 7(1 + 023 — k ) 



a. 



with 



a = f3(a + e)/(l +af3). On the other hand, if the core 
radius R is known, the wavenumber k of emitted waves 
can be analytically found if the condition Afc <C 1 is sat- 
isfied (cf. ||). The inverted dependence R(k) and the 
velocity V(k) are displayed in Fig. |l](a). 

For stationary pacemakers, the front velocity V van- 
ishes. This determines the wavenumber fco of a stationary 
pacemaker and thus allows us to find its core radius Rq. 
Analytical solutions for these key properties have been 
constructed. Wc find that 



fc = yl - cr/7 + az 3 ■ 
l + a/3 
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(3) 
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where fc„ lax = (a + e) {z\ — Z3 )/(/? — a) . The frequency flo 
of a stationary pacemaker is ilo = + (f3 — a)kQ. The 
wavenumber fco and the radius Rq of a stationary pace- 
maker are shown as functions of the coupling coefficient 
7 in Fig. |(b). 

Examining the constructed solutions, we note that gen- 
erally 1^3 < fio < The frequency Qq of a stationary 
pacemaker approaches the frequency Q 1 of rapid uni- 
form oscillations, when the core radius Rq — ► 00 (and 
^0 ~~ * k max ). On the other hand, when the core is small, 
fco is small and the frequency O is close to . Stationary 
pacemakers exist inside an interval of the coupling inten- 
sity 7 [see Fig. 0(b)] . Our approximate analysis based on 
singular perturbation theory is only valid when the core 
is not too small, i.e. Rq ;§> A. 

Some conclusions about the stability of stationary 
pacemakers can already be drawn from Fig. 0(a). Sup- 
pose the radius R has increased above the stationary ra- 
dius Rq. This leads to an increase of the wavenumber k of 
emitted waves which, in turn, will make the front velocity 
V negative. Therefore the front will retreat, decreasing 
the radius R back to its stationary value. This argument 
is, however, only applicable when the characteristic time 
scale of the core evolution is much longer than the time 
needed for the wave pattern to adjust to its changes, i.e. 
when t»1. Generally, the stability of stationary pace- 
makers should be numerically investigated. 

The system described by Eqs. (0) was integrated with 
an explicit Euler scheme where the Laplacian operator 
was discretized with a nearest-neighbor approximation. 
No-flux boundary conditions were used. Figure 0(a) dis- 
plays the evolution of a stationary pacemaker from a 
small initial perturbation of the real mode z. In the first 
stage, the core grows with approximately constant speed. 
Later, the growth is terminated and a stationary object 
is formed. Figure 0(b) displays the creation and emis- 
sion of waves in the oscillatory subsystem. The profile of 
the asymptotic stable stationary pacemaker is shown in 
Fig- |(c). 

Pacemakers are stable for sufficiently large r. When 
r is decreased, numerical integrations show that station- 
ary pacemakers become unstable. Close to the instabil- 
ity boundary, stable breathing and swinging pacemakers 
were found [Figs. 0(a)-0(b)] . For a breathing pacemaker, 
the center remains stationary whereas the radius oscil- 
lates. For swinging pacemakers, the radius stays approx- 
imately constant while the position of the pacemaker os- 
cillates. Further lowering of r leads to the disappearance 
of any stable pacemaker solutions. 

Similar to plane waves in the CGLE p9[ , the waves 
emitted by a pacemaker may become unstable when their 
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wavenumber kg exceeds the Eckhaus wavenumber given 
by fc E H ~ + z 3 )(l + a(3){3 + a/3 + 2a 2 )- 1 . Among 
other effects, this may iead to the destabilization of a 
stationary pacemaker, as illustrated by Figs. ||(c)-||(d). 
Phase singularities and thus oscillation amplitude defects 
are periodically generated at the core boundary, giving 
rise to a short-wavelength regime in the core and to a 
long-wavelength regime in the periphery. The core grad- 
ually grows and eventually the whole medium is occupied 
by rapid uniform oscillations. 

In contrast to pacemakers which are created by local 
heterogeneities, autonomous pacemakers are not pinned 
and their location is determined only by the initial con- 
ditions. Moreover, such self-organized structures are 
able to move through the medium when spatial param- 
eter gradients are present. Suppose, for example, that 
the parameter 7 varies with a constant gradient k, i.e. 
7(x) =70 + k(x — xq), where 70 is the value of 7 in 
the center of the pacemaker. For sufficiently small gradi- 
ents (kRq <C 7o)> linear perturbation theory can be used. 
Its application (see [p8| ) allows us to determine analyti- 
cally the drift velocity Vd as Vb = kRo 9 7 V(7o, ho) where 
V(7, k) is the front velocity given by Eq. (||). 

The simulation displayed in Fig. [|(a) was initiated with 
a stable stationary pacemaker. After a constant gradi- 
ent in the parameter 7 was introduced, the pacemaker 
drifted through the medium in the direction of increas- 
ing 7. When the gradient was removed, the drift of the 
pacemaker terminated and a spatially shifted stationary 
pacemaker was recovered. The emission and propagation 
of waves persisted during the drift [Fig. f|(b)]. In addi- 
tion, the Doppler effect led to a small increase of k in the 
direction of motion. 

In this Letter, we have analytically constructed self- 
organized pacemaker solutions in the vicinity of a pitch- 
fork- Hopf bifurcation. Our numerical investigations have 
shown that such self-organized patterns are stable for a 
wide range of parameters. To create autonomous pace- 
makers, a sufficiently strong local perturbation should be 
applied to the state corresponding to stable uniform os- 
cillations. This is in contrast to []15|, where autonomous 
target patterns were found near a Hopf bifurcation with 
a finite wavenumber and thus uniform oscillations of the 
medium were absolutely unstable. Our approach is also 
different from the model j|(J[n] that was constructed to 
explain target pattern formation in electrohydrodynamic 
convection and which is based on a Hopf bifurcation of 
a cellular spatial structure. On the other hand, Ohta 
et al. |3 have investigated a two-component activator- 
inhibitor model with coexistence of excitable kinetics 
and stable uniform oscillations, and reported several dif- 
ferent kinds of autonomous wave sources. The subse- 
quent numerical studies p(J have, however, shown that 
while localized target patterns are stable, target patterns 
which extend over the whole medium are unstable in the 
model and slowly evolve into uniform oscillations. Sta- 



ble localized target patterns are also found in the quintic 
CGLE || . 

Since our analysis is based on general amplitude 
equations, the results presented here are valid for any 
reaction-diffusion system near a soft onset of birhythmic- 
ity with small-amplitude limit cycles. In a separate pub- 
lication pj| , this analysis will be applied to a particular 
chemical model system. As in the case of a Turing-Hopf 
bifurcation, the results of our analysis based on the am- 
plitude equations may remain (qualitatively) applicable 
even at significant separation from the bifurcation point. 
Finally, we note that the physical mechanism responsi- 
ble for the stabilization of pacemakers in the considered 
system involves a long-range negative feedback, similar 
to the one necessary for the formation of stable local- 
ized spots in reaction-diffusion models with fast inhibitor 
diffusion. Here, however, an infinite-range inhibition is 
caused not by diffusion, but by non-damped propaga- 
tion of waves emitted from the core region. The effect of 
pacemaker drift in systems with spatial parameter gradi- 
ents provides a convenient experimental method to iden- 
tify self-organized pacemakers and distinguish them from 
other target patterns caused by local heterogeneities in 
the medium. 
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tion. Financial support of the Humboldt Foundation 
(Germany) is gratefully acknowledged. 
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FIG. 1. (a) Dependence of the front velocity V on the 
wavenumber k of emitted waves (solid line) and dependence 
of k on the core radius R (dashed line, plotted as R vs. k). 
(b) The wavenumber ko of the waves emitted by a station- 
ary pacemaker (solid line) and the corresponding radius Rq 
(dashed line) as functions of the coupling coefficient 7. The 
parameters are a — 1.4, /3 — 2.3, e — —2.1, A = 1, r = 5, 
v = 20, 7 = 0.13, cr = 0.1. 
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FIG. 2. Development of a stable stationary pacemaker 
(a,b) and its asymptotic profile (c). Frame (a) shows the 
evolution of the real mode amplitude z after an initial per- 
turbation. Frame (b) displays the evolution of ReA In frame 
(c) the spatial distribution of the variables z (solid line) , He A 
(dotted line) and \A\ (dashed line) are presented. The system 
size is L = 100 and the time interval is < t < 500; the same 
parameters as in Fig. [j]. In our gray-scale plots the black 
and white levels always correspond to the minimum and the 
maximum values of the plotted variable, respectively. In the 
space-time plots, time runs along the horizontal axis. 




FIG. 3. Breathing (a), swinging (b), and Eckhaus-unstable 
(c-d) pacemakers. The displayed coordinate and time ranges 
are AL = 50, AT = 125 (a-b) and AL = 100, AT = 625 
(c-d). The parameters are L — 100, a — 1.38, e = —3.18, 
A = 0.8, v = 83, 7 = 5.59 • 10" 4 , a = 3.4 ■ 10" 4 , and for 
(a): P = 3.0, t = 0.001, (b): f3 = 2.65, r = 0.001, (c-d): 
P = 2.1, r = 0.025. 
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FIG. 4. Drift of a pacemaker. The spatial gradient of the 
parameter 7 with /t/70 = 0.003 is applied inside the time 
interval indicated by vertical dashed lines in frame (a) , which 
shows the evolution of z in the time interval < t < 2 ■ 10 5 . 
The pacemaker is drifting in the direction of increased 7. The 
parameters are 70 = 5.59 • 10 - , (i = 2.3, r = 2. The rest are 
the same as in Fig. ^. Frame (b) displays the drifting wave 
pattern within a narrow time interval AT = 500 during the 
drift, marked by the dotted vertical line in frame (a). 
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